Visualizing Water: Interactive Timelapse Map of River Basins (draft)

draft imwi hydrology

Part 3: how to represent water flux over 3D maps of river basins.

Melanie BACOU https://linkedin/in/mbacou
2021-11-01

This notebook is Part 3 of an exploration to visualize results of hydrologic models. In Part 1 we built custom HTML widgets using D3.js, and in Part 2 we looked at rendering water fluxes using Sankey diagrams. Here we test multiple libraries to generate hillshade (3D) views of river basins and water infrastructure, in particular we want to compare leaflet, Three.js and MapboxGL implementations.

Aside from rendering topography and water streams in 3D (and potentially other covariate layers), our objective is to overlay custom labels to illustrate the hydrological cycle.

Another objective is to provide basin and sub-basin statistics, starting with precipitation, ET, soil moisture (used in the WA+ approach), and adding other covariates, such as population, land use allocation, and crop allocation.

Some inspiration below:

Show code
list.files("./fig", full.names=TRUE)[c(2,5)] %>%
  knitr::include_graphics()

WebGL: Sample Scene Processing with Rayshdaer

Data Acquisition

We’ll start with a sample scene of the Selingue Dam in the Niger River basin.

basin <- shapefile("~/Projects/2021-iwmi/data/mli/srtm/mli_basin.shp")
zoi <- shapefile("~/Projects/2021-iwmi/data/mli/srtm/zoi.shp")
ext <- extent(zoi)
center <- coordinates(zoi)

# Get CGIAR SRTM DEM at 90m
srtm <- getData("SRTM", lon=center[,1], lat=center[,2], path="_data") %>%
  crop(zoi) %>%
  mask(zoi)
srtm
class      : RasterLayer 
dimensions : 838, 973, 815374  (nrow, ncol, ncell)
resolution : 0.0008333333, 0.0008333333  (x, y)
extent     : -8.6725, -7.861667, 11.06917, 11.7675  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
source     : memory
names      : srtm_35_10 
values     : 326, 476  (min, max)
# Satellite basemaps
bmap <- maptiles::get_tiles(terra::ext(zoi), "Esri.WorldImagery", zoom=10) %>%
  stack() %>%
  crop(zoi) %>%
  mask(zoi)

The basin covers a large area, so we need 8 SRTM tiles, but 1 is enough for a proof of concept. Next we’ll get a satellite basemap.

plot(terra::vect(basin), col=pal[2], border=pal[1], lwd=2,
  main="Niger River Basin (Mali)")
plot(zoi, lty=3, col=alpha(pal["red"], .6), border=pal["red"], lwd=2, add=T)
text(-8, 10.5, "Selingue Dam\n(Mali)", col=pal["red"], cex=.7, font=2)
grid()

plot(terra::ext(zoi), 
  main="Selingue Dam (Mali) - ESRI World Imagery")
plotRGB(bmap, add=T)
grid(col="white")

plot(terra::rast(srtm),
  main="Selingue Dam (Mali) - SRTM 90m")
grid(col="white")

Scene Rendering

Next we convert the 2 rasters to a matrix format that’s compatible with Rayshader methods.

Show code
# Convert rasters to rayshader matrix format
srtm_array <- raster_to_matrix(srtm)

# Convert sat basemap to matrix (test)
r <- raster_to_matrix(bmap$red)
g <- raster_to_matrix(bmap$green)
b <- raster_to_matrix(bmap$blue)

bmap_array <- array(0, dim=c(nrow(r), ncol(r), 3))
bmap_array[,,1] <- r/255
bmap_array[,,2] <- g/255
bmap_array[,,3] <- b/255

bmap_array %>%
  aperm(c(2,1,3)) %>%
  # Stretch contrast
  rescale(to=c(0,1)) -> bmap_array
srtm_water <- srtm_array
srtm_water[srtm_water < 353] <- 0

basemap_sat <- srtm_array %>%
  height_shade() %>%
  add_overlay(bmap_array) %>%
  add_shadow(ray_shade(srtm_array, zscale=90)) %>%
  add_water(detect_water(srtm_water), color=alpha(pal["blue"], 0.4))

basemap_sat %>% plot_map()
Hillshade Basemap of Selingue Dam (Niger River basin)

Figure 1: Hillshade Basemap of Selingue Dam (Niger River basin)

That doesn’t look very clear, so instead we’ll create a basemap, not using the satellite image but a simple desert texture.

base_map <- srtm_array %>% 
  height_shade() %>% 
  add_overlay(sphere_shade(srtm_array, texture="desert", 
    zscale=4, colorintensity=5), alphalayer=0.5) %>%
  add_shadow(lamb_shade(srtm_array, zscale=6), 0) %>%
  add_shadow(ambient_shade(srtm_array), 0) %>%
  add_shadow(texture_shade(srtm_array, detail=8/10, contrast=9, brightness=11), 0.1) %>%
  add_water(detect_water(srtm_water), color=alpha(pal["blue"], 0.4))

saveRDS(base_map, "./_data/base_map.rds")
Show code
# Upsampled SRTM
alt <- getData("alt", country="MLI", path="_data") %>%
  crop(zoi) %>%
  raster_to_matrix()

base_map <- readRDS("./_data/base_map.rds")

base_map %>% plot_3d(alt,
  zscale=20, zoom=0.5,
  phi=30, theta=20, fov=0,
  shadow=TRUE, shadowcolor=pal["black"])

rglwidget()

Figure 2: Step 1 - Hillshade Basemap of Selingue Dam

Looks better, so let’s acquire and overlay spatial features from OSM.

osm_roads <- opq(bbox(zoi)) %>% 
  add_osm_feature("highway") %>% 
  osmdata_sf()

osm_water = opq(bbox(zoi)) %>% 
  add_osm_feature("waterway") %>% 
  osmdata_sf()

osm_place = opq(bbox(zoi)) %>% 
  add_osm_feature("place") %>% 
  osmdata_sf()

road_layer <- generate_line_overlay(
  dplyr::filter(osm_roads$osm_lines, highway %in% c("primary", "secondary")),
  extent=ext, srtm_array, linewidth=5, color=pal["black"])

water_layer <- generate_line_overlay(
  osm_water$osm_lines, 
  extent=ext, srtm_array, linewidth=3, color=pal["blue"])

place_layer <- generate_label_overlay(
  dplyr::filter(osm_place$osm_points, !is.na(name) & nchar(name)<10), 
  extent=ext, heightmap=srtm_array,
  font=2, text_size=1.6, point_size=1.6, color=pal["black"],
  halo_color="white", halo_expand=2, halo_blur=1, halo_alpha=.9, seed=1,
  data_label_column="name")

scene <- base_map %>% 
#scene <- basemap_sat %>%   
  add_overlay(road_layer) %>%
  add_overlay(water_layer, alphalayer=1) %>%
  add_overlay(place_layer)

saveRDS(scene, "./_data/scene.rds")

Finally we’ll use WebGL to render this scene in 3D.

amb_layer <- ambient_shade(srtm_array, zscale=1/5)

scene2 <- srtm_array %>% 
  height_shade() %>%
  add_shadow(texture_shade(srtm_array, detail=8/10, contrast=9, brightness=11), 0) %>%
  add_shadow(amb_layer, 0) %>%
  add_overlay(road_layer) %>%
  add_overlay(water_layer, alphalayer=1) %>%
  add_overlay(place_layer)

saveRDS(scene2, "./_data/scene2.rds")
Show code
scene <- readRDS("./_data/scene.rds")
scene2 <- readRDS("./_data/scene2.rds")

scene %>% plot_3d(alt, zscale=20, 
  theta=30, phi=20, fov=0, zoom=0.5)

# Add polygon annotations
xyz <- sf::read_sf("~/Projects/2021-iwmi/data/mli/srtm/xyz.shp")
render_polygons(xyz, ext, data_column_top="z",
  scale_data=1, color=alpha(pal["orange"], 0.8), 
  lit=F, light_intensity=0.01, clear_previous=T)

rglwidget()

Figure 3: Step 2: Interactive 3D Scene of Selingue Dam

WebGL: Full Scene with Rayshader

Following the approach above, we generate a 3D scene for the entire Niger River basin used in the analysis.

Data Acquisition

Show code
# Basin boundaries
zoi <- shapefile("~/Projects/2021-iwmi/data/mli/srtm/mli_basin.shp")
ext <- extent(zoi)
center <- coordinates(zoi)

# Admin2 boundaries
adm <- sf::st_read(
  "~/Projects/2021-iwmi/data/mli/srtm/mli_adm2_lines.shp")
places <- sf::st_read(
  "~/Projects/2021-iwmi/data/mli/srtm/mli_adm2_centroids.shp")

# Get upsampled CGIAR SRTM DEM at 90m
alt <- lapply(c("MLI", "GIN", "CIV"), 
  function(x) getData("alt", country=x, mask=FALSE) %>%
    crop(zoi) %>%
    mask(zoi)
)
alt <- mosaic(alt[[1]], alt[[2]], alt[[3]], fun=mean)
alt <- mask(crop(alt, zoi), zoi)
writeRaster(alt, "~/Projects/2021-iwmi/data/mli/srtm/mli_alt.tif",
  overwrite=T)
alt

# Get ESA WorldCover (clipped in QGIS)
luc <- raster("~/Projects/2021-iwmi/data/mli/srtm/ESA WorldCover clip.tiff")

# Note that I can't find a colormap for this raster, so switching so ESA CCI
luc <- raster("~/Projects/2021-iwmi/data/mli/srtm/mli_esa_cci_300m.tif")

# And its colormap
pal.luc <- fread("~/Maps/ESA/ESA CCI Colormap.txt")
setnames(pal.luc, c("value", "R", "G", "B", "A", "label"))
pal.luc[, hex := rgb(R, G, B, maxColorValue=255)]

# OSM basemap
bmap <- maptiles::get_tiles(terra::ext(zoi), "OpenStreetMap.HOT", zoom=6) %>%
  stack() %>%
  crop(zoi) %>%
  mask(zoi)

# We'll use GloRIC v10 stream network instead of OSM
# Features are filtered to 'Class_hydr < 12'
gloric <- sf::st_read("~/Projects/2021-iwmi/data/mli/srtm/mli_gloric_filtered.shp")

save(zoi, ext, adm, places, alt, luc, pal.luc, bmap, gloric,
  file="./_data/osm.RData")
load("./_data/osm.RData")

plot(terra::rast(alt),
  main="Niger River Basin (Mali) - SRTM 1km")
plot(adm, lty=3, col=pal["black"], lwd=1, add=T)
plot(gloric, col=pal["blue"], lwd=.4, add=T)
grid()

plot(terra::rast(luc), col=pal.luc$hex, 
  breaks=pal.luc$value, legend=F,
  main="Niger River Basin (Mali) - ESA CCI 300m")
plot(adm, lty=3, col=pal["black"], lwd=1, add=T)
grid()

Scene Rendering

Same as above.

Show code
# Convert rasters to rayshader matrix format
srtm_array <- raster_to_matrix(alt)

# Convert sat basemap to matrix (test)
r <- raster_to_matrix(bmap$red)
g <- raster_to_matrix(bmap$green)
b <- raster_to_matrix(bmap$blue)

bmap_array <- array(0, dim=c(nrow(r), ncol(r), 3))
bmap_array[,,1] <- r/255
bmap_array[,,2] <- g/255
bmap_array[,,3] <- b/255

bmap_array %>%
  aperm(c(2,1,3)) %>%
  # Stretch contrast
  rescale(to=c(0,1)) -> bmap_array
base_map <- srtm_array %>% 
  height_shade() %>% 
  add_overlay(sphere_shade(srtm_array, texture="desert"), alphalayer=0.5) %>%
  add_shadow(lamb_shade(srtm_array, zscale=6), 0) %>%
  add_shadow(ambient_shade(srtm_array), 0)

saveRDS(base_map, "./_data/base_map_basin.rds")
Show code
base_map <- readRDS("./_data/base_map_basin.rds")
plot_map(base_map)
Step 1 - Hillshade Basemap Niger River Basin

Figure 4: Step 1 - Hillshade Basemap Niger River Basin

admin_layer <- generate_line_overlay(adm,
  extent=ext, srtm_array, linewidth=3, lty=1, color=pal["light"])

water_layer <- generate_line_overlay(gloric, 
  extent=ext, srtm_array, linewidth=3, color=pal["light-blue"])

place_layer <- generate_label_overlay(places,
  extent=ext, heightmap=srtm_array,
  font=2, text_size=1.6, point_size=0, color=pal["light"],
  halo_color=pal["black"], halo_expand=2, halo_blur=2, halo_alpha=.6, seed=1,
  data_label_column="ADM2_NAME")

scene <- base_map %>% 
#scene <- basemap_luc %>%   
  add_overlay(admin_layer) %>%
  add_overlay(water_layer)
  #add_overlay(place_layer)

saveRDS(scene, "./_data/scene_basin.rds")

Finally we’ll use WebGL to render this scene in 3D.

Show code
try(rgl.close())
scene <- readRDS("./_data/scene_basin.rds")

scene %>% plot_3d(resize_matrix(srtm_array, 1/3), zscale=80, 
  theta=30, phi=50, fov=0, zoom=0.5)

rglwidget()

Figure 5: Step 2: Interactive 3D Scene of Niger River Basin

3D Rendering with Three.js

Similar implementation using Three.js libraries and viewer. The scene is 12 MB in size (low resampling), it was generated using the Qgis2threejs plugin. The main difference is we’re showing ESA CCI land cover classification as an image overlay.

Akagi, Minoru. 2021. Qgis2threejs Plugin 2.4 documentation.” https://qgis2threejs.readthedocs.io/en/docs/index.html.
Franke, Loraine, and Daniel Haehn. 2020. Modern scientific visualizations on the web.” Informatics 7 (4). https://doi.org/10.3390/INFORMATICS7040037.
Hssaisoune, Mohammed, Lhoussaine Bouchaou, Abdelfattah Sifeddine, Ilham Bouimetarhan, and Abdelghani Chehbouni. 2020. Moroccan Groundwater Resources and Evolution with Global Climate Changes.” Geosciences (Switzerland) 10 (2). https://doi.org/10.3390/geosciences10020081.
Liu, Li, Deborah Silver, and Karen Bemis. 2019. Visualizing three-dimensional ocean eddies in web browsers.” IEEE Access 7: 44734–47. https://doi.org/10.1109/ACCESS.2019.2909655.
Marre, Alain. 2015. Eau, Agriculture et Pauvreté dans le bassin du Niger: synthèse des résultats du BFP Niger.” Physio-Géo, no. Volume 9 (January): 21–24. https://doi.org/10.4000/physio-geo.4650.
Morgan, Tyler. 2021. A Step-by-Step Guide to Making 3D Maps with Satellite Imagery in R.” https://www.tylermw.com/a-step-by-step-guide-to-making-3d-maps-with-satellite-imagery-in-r/.
Popovs, Konrāds, Tomas Saks, and Jānis Jātnieks. 2015. Web Based 3D Visualization of Kazu Leja Conceptual Geological Model.” Estonian Journal of Earth Sciences 64 (2): 173–88. https://doi.org/10.3176/earth.2015.25.

References

Corrections

If you see mistakes or want to suggest changes, please create an issue on the source repository.

Reuse

Text and figures are licensed under Creative Commons Attribution CC BY 4.0. Source code is available at https://github.com/mbacou/mb-labs, unless otherwise noted. The figures that have been reused from other sources don't fall under this license and can be recognized by a note in their caption: "Figure from ...".

Citation

For attribution, please cite this work as

BACOU (2021, Nov. 1). Mel's Labs: Visualizing Water: Interactive Timelapse Map of River Basins (draft). Retrieved from https://mbacou.github.io/mb-labs/posts/2021-10-24-3dmap/

BibTeX citation

@misc{bacou2021visualizing,
  author = {BACOU, Melanie},
  title = {Mel's Labs: Visualizing Water: Interactive Timelapse Map of River Basins (draft)},
  url = {https://mbacou.github.io/mb-labs/posts/2021-10-24-3dmap/},
  year = {2021}
}